/*******************************************************************************
 * Copyright 2009, 2017 Martin Davis
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

package org.locationtech.proj4j.util;

import org.locationtech.proj4j.ConvergenceFailureException;
import org.locationtech.proj4j.InvalidValueException;


public class ProjectionMath {

    public final static double PI = Math.PI;

    public final static double HALFPI = Math.PI / 2.0;

    public final static double QUARTERPI = Math.PI / 4.0;

    public final static double FORTPI = QUARTERPI;

    public final static double TWOPI = Math.PI * 2.0;

    public final static double RTD = 180.0 / Math.PI;
    
    public final static double DTR = Math.PI / 180.0;

    public final static double EPS10 = 1.0e-10;

    /**
     * Degree versions of trigonometric functions
     */
    public static double sind(double v) {
        return Math.sin(v * DTR);
    }

    public static double cosd(double v) {
        return Math.cos(v * DTR);
    }

    public static double tand(double v) {
        return Math.tan(v * DTR);
    }

    public static double asind(double v) {
        return Math.asin(v) * RTD;
    }

    public static double acosd(double v) {
        return Math.acos(v) * RTD;
    }

    public static double atand(double v) {
        return Math.atan(v) * RTD;
    }

    public static double atan2d(double y, double x) {
        return Math.atan2(y, x) * RTD;
    }

    public static double asin(double v) {
        if (Math.abs(v) > 1.)
            return v < 0.0 ? -Math.PI / 2 : Math.PI / 2;
        return Math.asin(v);
    }

    public static double acos(double v) {
        if (Math.abs(v) > 1.)
            return v < 0.0 ? Math.PI : 0.0;
        return Math.acos(v);
    }

    public static double sqrt(double v) {
        return v < 0.0 ? 0.0 : Math.sqrt(v);
    }

    public static double distance(double dx, double dy) {
        return Math.sqrt(dx * dx + dy * dy);
    }

    public static double hypot(double x, double y) {
        if (x < 0.0)
            x = -x;
        else if (x == 0.0)
            return y < 0.0 ? -y : y;
        if (y < 0.0)
            y = -y;
        else if (y == 0.0)
            return x;
        if (x < y) {
            x /= y;
            return y * Math.sqrt(1.0 + x * x);
        } else {
            y /= x;
            return x * Math.sqrt(1.0 + y * y);
        }
    }

    public static double atan2(double y, double x) {
        return Math.atan2(y, x);
    }

    public static double trunc(double v) {
        return v < 0.0 ? Math.ceil(v) : Math.floor(v);
    }

    public static double frac(double v) {
        return v - trunc(v);
    }

    public static double degToRad(double v) {
        return v * Math.PI / 180.0;
    }

    public static double radToDeg(double v) {
        return v * 180.0 / Math.PI;
    }

    // For negative angles, d should be negative, m & s positive.
    public static double dmsToRad(double d, double m, double s) {
        if (d >= 0)
            return (d + m / 60 + s / 3600) * Math.PI / 180.0;
        return (d - m / 60 - s / 3600) * Math.PI / 180.0;
    }

    // For negative angles, d should be negative, m & s positive.
    public static double dmsToDeg(double d, double m, double s) {
        if (d >= 0)
            return (d + m / 60 + s / 3600);
        return (d - m / 60 - s / 3600);
    }

    public static double normalizeLatitude(double angle) {
        if (Double.isInfinite(angle) || Double.isNaN(angle))
            throw new InvalidValueException("Infinite latitude");
        while (angle > ProjectionMath.HALFPI)
            angle -= Math.PI;
        while (angle < -ProjectionMath.HALFPI)
            angle += Math.PI;
        return angle;
//		return Math.IEEEremainder(angle, Math.PI);
    }

    public static double normalizeLongitude(double angle) {
        if (Double.isInfinite(angle) || Double.isNaN(angle))
            throw new InvalidValueException("Infinite longitude");
        while (angle > Math.PI)
            angle -= TWOPI;
        while (angle < -Math.PI)
            angle += TWOPI;
        return angle;
//		return Math.IEEEremainder(angle, Math.PI);
    }

    public static double normalizeAngle(double angle) {
        if (Double.isInfinite(angle) || Double.isNaN(angle))
            throw new InvalidValueException("Infinite angle");
        while (angle > TWOPI)
            angle -= TWOPI;
        while (angle < 0)
            angle += TWOPI;
        return angle;
    }

/*
	public static void latLongToXYZ(Point2D.Double ll, Point3D xyz) {
		double c = Math.cos(ll.y);
		xyz.x = c * Math.cos(ll.x);
		xyz.y = c * Math.sin(ll.x);
		xyz.z = Math.sin(ll.y);
	}

	public static void xyzToLatLong(Point3D xyz, Point2D.Double ll) {
		ll.y = MapMath.asin(xyz.z);
		ll.x = MapMath.atan2(xyz.y, xyz.x);
	}
*/

    public static double greatCircleDistance(double lon1, double lat1, double lon2, double lat2) {
        double dlat = Math.sin((lat2 - lat1) / 2);
        double dlon = Math.sin((lon2 - lon1) / 2);
        double r = Math.sqrt(dlat * dlat + Math.cos(lat1) * Math.cos(lat2) * dlon * dlon);
        return 2.0 * Math.asin(r);
    }

    public static double sphericalAzimuth(double lat0, double lon0, double lat, double lon) {
        double diff = lon - lon0;
        double coslat = Math.cos(lat);

        return Math.atan2(
                coslat * Math.sin(diff),
                (Math.cos(lat0) * Math.sin(lat) -
                        Math.sin(lat0) * coslat * Math.cos(diff))
        );
    }

    public static boolean sameSigns(double a, double b) {
        return a < 0 == b < 0;
    }

    public static boolean sameSigns(int a, int b) {
        return a < 0 == b < 0;
    }

    public static double takeSign(double a, double b) {
        a = Math.abs(a);
        if (b < 0)
            return -a;
        return a;
    }

    public static int takeSign(int a, int b) {
        a = Math.abs(a);
        if (b < 0)
            return -a;
        return a;
    }
/*
  public static double distance(Point2D.Double a, Point2D.Double b) {
    return distance(a.x-b.x, a.y-b.y);
  }

	public final static int DONT_INTERSECT = 0;
	public final static int DO_INTERSECT = 1;
	public final static int COLLINEAR = 2;

	public static int intersectSegments(Point2D.Double aStart, Point2D.Double aEnd, Point2D.Double bStart, Point2D.Double bEnd, Point2D.Double p) {
		double a1, a2, b1, b2, c1, c2;
		double r1, r2, r3, r4;
		double denom, offset, num;

		a1 = aEnd.y-aStart.y;
		b1 = aStart.x-aEnd.x;
		c1 = aEnd.x*aStart.y - aStart.x*aEnd.y;
		r3 = a1*bStart.x + b1*bStart.y + c1;
		r4 = a1*bEnd.x + b1*bEnd.y + c1;

		if (r3 != 0 && r4 != 0 && sameSigns(r3, r4))
			return DONT_INTERSECT;

		a2 = bEnd.y-bStart.y;
		b2 = bStart.x-bEnd.x;
		c2 = bEnd.x*bStart.y-bStart.x*bEnd.y;
		r1 = a2*aStart.x + b2*aStart.y + c2;
		r2 = a2*aEnd.x + b2*aEnd.y + c2;

		if (r1 != 0 && r2 != 0 && sameSigns(r1, r2))
			return DONT_INTERSECT;

		denom = a1*b2 - a2*b1;
		if (denom == 0)
			return COLLINEAR;

		offset = denom < 0 ? -denom/2 : denom/2;

		num = b1*c2 - b2*c1;
		p.x = (num < 0 ? num-offset : num+offset) / denom;

		num = a2*c1 - a1*c2;
		p.y = (num < 0 ? num-offset : num+offset) / denom;

		return DO_INTERSECT;
	}

  /*
	public static double dot(Point2D.Double a, Point2D.Double b) {
		return a.x*b.x + a.y*b.y;
	}
	
	public static Point2D.Double perpendicular(Point2D.Double a) {
		return new Point2D.Double(-a.y, a.x);
	}
	
	public static Point2D.Double add(Point2D.Double a, Point2D.Double b) {
		return new Point2D.Double(a.x+b.x, a.y+b.y);
	}
	
	public static Point2D.Double subtract(Point2D.Double a, Point2D.Double b) {
		return new Point2D.Double(a.x-b.x, a.y-b.y);
	}
	
	public static Point2D.Double multiply(Point2D.Double a, Point2D.Double b) {
		return new Point2D.Double(a.x*b.x, a.y*b.y);
	}
	
	public static double cross(Point2D.Double a, Point2D.Double b) {
		return a.x*b.y - b.x*a.y;
	}
  
  public static void normalize(Point2D.Double a) {
    double d = distance(a.x, a.y);
    a.x /= d;
    a.y /= d;
  }
  
  public static void negate(Point2D.Double a) {
    a.x = -a.x;
    a.y = -a.y;
  }
  

*/

    public static double cross(double x1, double y1, double x2, double y2) {
        return x1 * y2 - x2 * y1;
    }

    public static double longitudeDistance(double l1, double l2) {
        return Math.min(
                Math.abs(l1 - l2),
                ((l1 < 0) ? l1 + Math.PI : Math.PI - l1) + ((l2 < 0) ? l2 + Math.PI : Math.PI - l2)
        );
    }

    public static double geocentricLatitude(double lat, double flatness) {
        double f = 1.0 - flatness;
        return Math.atan((f * f) * Math.tan(lat));
    }

    public static double geographicLatitude(double lat, double flatness) {
        double f = 1.0 - flatness;
        return Math.atan(Math.tan(lat) / (f * f));
    }

    public static double tsfn(double phi, double sinphi, double e) {
        sinphi *= e;
        return (Math.tan(.5 * (ProjectionMath.HALFPI - phi)) /
                Math.pow((1. - sinphi) / (1. + sinphi), .5 * e));
    }

    public static double msfn(double sinphi, double cosphi, double es) {
        return cosphi / Math.sqrt(1.0 - es * sinphi * sinphi);
    }

    private final static int N_ITER = 15;

    public static double phi2(double ts, double e) {
        double eccnth, phi, con, dphi;
        int i;

        eccnth = .5 * e;
        phi = ProjectionMath.HALFPI - 2. * Math.atan(ts);
        i = N_ITER;
        do {
            con = e * Math.sin(phi);
            dphi = ProjectionMath.HALFPI - 2. * Math.atan(ts * Math.pow((1. - con) / (1. + con), eccnth)) - phi;
            phi += dphi;
        } while (Math.abs(dphi) > 1e-10 && --i != 0);
        if (i <= 0)
            throw new ConvergenceFailureException("Computation of phi2 failed to converage after " + N_ITER + " iterations");
        return phi;
    }

    private final static double C00 = 1.0;
    private final static double C02 = .25;
    private final static double C04 = .046875;
    private final static double C06 = .01953125;
    private final static double C08 = .01068115234375;
    private final static double C22 = .75;
    private final static double C44 = .46875;
    private final static double C46 = .01302083333333333333;
    private final static double C48 = .00712076822916666666;
    private final static double C66 = .36458333333333333333;
    private final static double C68 = .00569661458333333333;
    private final static double C88 = .3076171875;
    private final static int MAX_ITER = 10;

    public static double[] enfn(double es) {
        double t;
        double[] en = new double[5];
        en[0] = C00 - es * (C02 + es * (C04 + es * (C06 + es * C08)));
        en[1] = es * (C22 - es * (C04 + es * (C06 + es * C08)));
        en[2] = (t = es * es) * (C44 - es * (C46 + es * C48));
        en[3] = (t *= es) * (C66 - es * C68);
        en[4] = t * es * C88;
        return en;
    }

    public static double mlfn(double phi, double sphi, double cphi, double[] en) {
        cphi *= sphi;
        sphi *= sphi;
        return en[0] * phi - cphi * (en[1] + sphi * (en[2] + sphi * (en[3] + sphi * en[4])));
    }

    public static double inv_mlfn(double arg, double es, double[] en) {
        double s, t, phi, k = 1. / (1. - es);

        phi = arg;
        for (int i = MAX_ITER; i != 0; i--) {
            s = Math.sin(phi);
            t = 1. - es * s * s;
            phi -= t = (mlfn(phi, s, Math.cos(phi), en) - arg) * (t * Math.sqrt(t)) * k;
            if (Math.abs(t) < 1e-11)
                return phi;
        }
        return phi;
    }

    private final static double P00 = .33333333333333333333;
    private final static double P01 = .17222222222222222222;
    private final static double P02 = .10257936507936507936;
    private final static double P10 = .06388888888888888888;
    private final static double P11 = .06640211640211640211;
    private final static double P20 = .01641501294219154443;

    public static double[] authset(double es) {
        double t;
        double[] APA = new double[3];
        APA[0] = es * P00;
        t = es * es;
        APA[0] += t * P01;
        APA[1] = t * P10;
        t *= es;
        APA[0] += t * P02;
        APA[1] += t * P11;
        APA[2] = t * P20;
        return APA;
    }

    public static double authlat(double beta, double[] APA) {
        double t = beta + beta;
        return (beta + APA[0] * Math.sin(t) + APA[1] * Math.sin(t + t) + APA[2] * Math.sin(t + t + t));
    }

    public static double qsfn(double sinphi, double e, double one_es) {
        double con;

        if (e >= 1.0e-7) {
            con = e * sinphi;
            return (one_es * (sinphi / (1. - con * con) -
                    (.5 / e) * Math.log((1. - con) / (1. + con))));
        } else
            return (sinphi + sinphi);
    }

    /*
     * Java translation of "Nice Numbers for Graph Labels"
     * by Paul Heckbert
     * from "Graphics Gems", Academic Press, 1990
     */
    public static double niceNumber(double x, boolean round) {
        int expv;				/* exponent of x */
        double f;				/* fractional part of x */
        double nf;				/* nice, rounded fraction */

        expv = (int) Math.floor(Math.log(x) / Math.log(10));
        f = x / Math.pow(10., expv);		/* between 1 and 10 */
        if (round) {
            if (f < 1.5)
                nf = 1.;
            else if (f < 3.)
                nf = 2.;
            else if (f < 7.)
                nf = 5.;
            else
                nf = 10.;
        } else if (f <= 1.)
            nf = 1.;
        else if (f <= 2.)
            nf = 2.;
        else if (f <= 5.)
            nf = 5.;
        else
            nf = 10.;
        return nf * Math.pow(10., expv);
    }

    /**
     * Evaluate complex polynomial.
     * Note coefficients are always C[1] to C[n], C[0] is always (0,0).
     */
    public static Complex zpoly1(Complex z, Complex[] c) {
        Complex a = new Complex(c[c.length - 1]);
        double t;
        int n = c.length - 1;
        while (n-- > 0) {
            Complex C = c[n];
            a.r = C.r + z.r * (t = a.r) - z.i * a.i;
            a.i = C.i + z.r * a.i + z.i * t;
        }
        a.r = z.r * (t = a.r) - z.i * a.i;
        a.i = z.r * a.i + z.i * t;
        return a;
    }

    /**
     * Evaluate a complex polynomial and its derivative
     */
    public static Complex zpoly1d(Complex z, Complex[] C, Complex der) {
        Complex a, b;
        double t;
        boolean first = true;

        a = new Complex(C[C.length - 1]);
        b = new Complex(a);
        for (int i = C.length - 1; i > 0; i--) {
            if (first) {
                first = false;
            } else {
                b.r = a.r + z.r * (t = b.r) - z.i * b.i;
                b.i = a.i + z.r * b.i + z.i * t;
            }
            Complex c = C[i];
            a.r = c.r + z.r * (t = a.r) - z.i * a.i;
            a.i = c.i + z.r * a.i + z.i * t;
        }

        b.r = a.r + z.r * (t = b.r) - z.i * b.i;
        b.i = a.i + z.r * b.i + z.i * t;
        a.r = z.r * (t = a.r) - z.i * a.i;
        a.i = z.r * a.i + z.i * t;
        der.i = b.i;
        der.r = b.r;
        return a;
    }

    /**
     * Tests whether the datum parameter-based transform
     * is the identity transform
     * (in which case datum transformation can be short-circuited,
     * thus avoiding some loss of numerical precision).
     *
     * @param transform
     * @return
     */
    public static boolean isIdentity(double[] transform) {
        for (int i = 0; i < transform.length; i++) {
            // scale factor will normally be 1 for an identity transform
            if (i == 6) {
                if (transform[i] != 1.0 && transform[i] != 0.0)
                    return false;
            } else if (transform[i] != 0.0) return false;
        }
        return true;
    }


    /* SECONDS_TO_RAD = Pi/180/3600 */
    public static final double SECONDS_TO_RAD = 4.84813681109535993589914102357e-6;
    public static final double MILLION = 1000000.0;
}
